home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / os2 / sysb091a.zip / sysbench / src / pmb_fft.c < prev    next >
C/C++ Source or Header  |  1994-11-05  |  12KB  |  573 lines
  1. /*
  2.  C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
  3.  by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
  4. */
  5.  
  6. //#include <stdio.h>
  7. #include <math.h>
  8.  
  9. /***************************************************************/
  10. /* Timer options. You MUST uncomment one of the options below  */
  11. /* or compile, for example, with the '-DUNIX' option.          */
  12. /*                                                             */
  13. /* Example compilation:                                        */
  14. /* UNIX systems:                                               */
  15. /* cc -DUNIX -O tfftdp.c -o fft -lm                            */
  16. /***************************************************************/
  17. /* #define Amiga       */
  18. /* #define UNIX        */
  19. /* #define UNIX_Old    */
  20. /* #define VMS         */
  21. /* #define BORLAND_C   */
  22. /* #define MSC         */
  23. /* #define MAC         */
  24. /* #define IPSC        */
  25. /* #define FORTRAN_SEC */
  26. /* #define GTODay      */
  27. /* #define CTimer      */
  28. /* #define UXPM        */
  29.  
  30. #ifndef PI
  31. #define PI  3.14159265358979323846
  32. #endif
  33.  
  34. static double xr[262144],xi[262144];
  35.  
  36. double pmb_fft()
  37. {
  38. int i,j,np,npm,lp,n2,kr,ki;
  39. double *pxr,*pxi;
  40. double a,enp,t,x,y,z,zr,zi,zm;
  41. double t1,t2,t3,benchtime,VAX_REF,dtime();
  42. int dfft();
  43.  
  44. np  = 8;
  45. pxr = xr;
  46. pxi = xi;
  47.  
  48. //      printf("\nFFT benchmark - Double Precision - V1.0 -  05 Jan 1993\n\n");
  49.  
  50. //      printf("FFT size   Time(sec)   max error\n");
  51.  
  52.       t3 = dtime();
  53.       for (lp = 1; lp <= 15; lp++ )
  54.       { 
  55.       np = np * 2; 
  56. //      printf(" %7d ",np);
  57.       enp = np; 
  58.       npm = np / 2  - 1;  
  59.       t = PI / enp;
  60.       *pxr = (enp - 1.0) * 0.5;
  61.       *pxi = 0.0;
  62.       n2 = np / 2;  
  63.       *(pxr+n2) = -0.5;
  64.       *(pxi+n2) =  0.0;
  65.       
  66.       for (i = 1; i <= npm; i++) 
  67.       {
  68.           j = np - i;
  69.           *(pxr+i) = -0.5;
  70.           *(pxr+j) = -0.5;
  71.           z = t * (double)i;  
  72.           y = -0.5*(cos(z)/sin(z));
  73.           *(pxi+i) =  y;
  74.           *(pxi+j) = -y;
  75.       }
  76.       
  77.       t1 = dtime();
  78.       dfft(xr,xi,np);
  79.       t2 = dtime() - t1;
  80.       if (t2 < 0.0) t2 = 0.0;
  81. //      printf(" %10.4lf ",t2);
  82. /*
  83.       for (i=0; i<=15; i++) printf("%d  %g  %g\n",i,xr[i],xi[i]);
  84. */
  85.       zr = 0.0; 
  86.       zi = 0.0; 
  87.       kr = 0; 
  88.       ki = 0; 
  89.       npm = np-1;
  90.       for (i = 0; i <= npm; i++ ) 
  91.       { 
  92.           a = fabs( xr[i] - i );
  93.           if (zr < a) 
  94.           { 
  95.          zr = a; 
  96.          kr = i; 
  97.           }
  98.           
  99.           a = fabs(xi[i]);
  100.           if (zi < a) 
  101.           { 
  102.          zi = a; 
  103.          ki = i; 
  104.           }
  105.       }
  106.       zm = zr;
  107.       if ( fabs(zr) < fabs(zi) ) zm = zi;
  108. //      printf(" %10.2g %10.4\n",zm,t2);
  109.  
  110. /*
  111.       printf("re %7d  %10.2g  im %7d  %10.2g\n",kr,zr,ki,zi);
  112. */
  113.       }
  114.       benchtime = dtime() - t3;
  115.       VAX_REF = 140.658;
  116. //      printf("\nBenchTime (sec) = %10.4lf\n",benchtime);
  117. //      printf("VAX_FFTs = %10.3lf\n\n",VAX_REF/benchtime);
  118.       return VAX_REF/benchtime;
  119. }
  120.  
  121.  
  122.  
  123. /********************************************************/
  124. /*     A Duhamel-Hollman split-radix dif fft            */
  125. /*     Ref: Electronics Letters, Jan. 5, 1984           */
  126. /*     Complex input and output data in arrays x and y  */
  127. /*     Length is n.                                     */
  128. /********************************************************/
  129.  
  130. static int dfft( x, y, np )
  131. double x[]; double y[]; int np ;
  132. {
  133. double *px,*py;
  134. int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4;
  135. double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt,tpi;
  136.  
  137.   px = x - 1; 
  138.   py = y - 1;
  139.   i = 2; 
  140.   m = 1; 
  141.   
  142.   while (i < np) 
  143.   { 
  144.   i = i+i; 
  145.   m = m+1; 
  146.   }
  147.   
  148.   n = i; 
  149.   
  150.   if (n != np) 
  151.   {
  152.      for (i = np+1; i <= n; i++)  
  153.      { 
  154.      *(px + i) = 0.0; 
  155.      *(py + i) = 0.0; 
  156.      }
  157.      printf("nuse %d point fft",n); 
  158.   }
  159.  
  160.   n2 = n+n;
  161.   tpi = 2.0 * PI;
  162.   for (k = 1;  k <= m-1; k++ ) 
  163.   {
  164.     n2 = n2 / 2; 
  165.     n4 = n2 / 4; 
  166.     e  = tpi / (double)n2; 
  167.     a  = 0.0;
  168.  
  169.     for (j = 1; j<= n4 ; j++) 
  170.     {
  171.       a3 = 3.0 * a; 
  172.       cc1 = cos(a); 
  173.       ss1 = sin(a);
  174.       
  175.       cc3 = cos(a3); 
  176.       ss3 = sin(a3); 
  177.       a = e * (double)j; 
  178.       is = j; 
  179.       id = 2 * n2;
  180.       
  181.       while ( is < n ) 
  182.       {
  183.          for (i0 = is; i0 <= n-1; i0 = i0 + id) 
  184.          {
  185.          i1 = i0 + n4; 
  186.          i2 = i1 + n4; 
  187.          i3 = i2 + n4;
  188.          r1 = *(px+i0) - *(px+i2);
  189.          *(px+i0) = *(px+i0) + *(px+i2);
  190.          r2 = *(px+i1) - *(px+i3);
  191.          *(px+i1) = *(px+i1) + *(px+i3);
  192.          s1 = *(py+i0) - *(py+i2);
  193.          *(py+i0) = *(py+i0) + *(py+i2);
  194.          s2 = *(py+i1) - *(py+i3);
  195.          *(py+i1) = *(py+i1) + *(py+i3);
  196.          s3 = r1 - s2; r1 = r1 + s2; 
  197.          s2 = r2 - s1; r2 = r2 + s1;
  198.          *(px+i2) = r1*cc1 - s2*ss1; 
  199.          *(py+i2) = -s2*cc1 - r1*ss1;
  200.          *(px+i3) = s3*cc3 + r2*ss3;
  201.          *(py+i3) = r2*cc3 - s3*ss3;
  202.          }
  203.       is = 2 * id - n2 + j; 
  204.       id = 4 * id;
  205.       }
  206.     }
  207.   }
  208.  
  209. /************************************/
  210. /*  Last stage, length=2 butterfly  */
  211. /************************************/
  212.   is = 1; 
  213.   id = 4;
  214.   
  215.   while ( is < n) 
  216.   {
  217.      for (i0 = is; i0 <= n; i0 = i0 + id) 
  218.      {
  219.       i1 = i0 + 1; 
  220.       r1 = *(px+i0);
  221.       *(px+i0) = r1 + *(px+i1);
  222.       *(px+i1) = r1 - *(px+i1);
  223.       r1 = *(py+i0);
  224.       *(py+i0) = r1 + *(py+i1);
  225.       *(py+i1) = r1 - *(py+i1);
  226.      }
  227.   is = 2*id - 1; 
  228.   id = 4 * id; 
  229.   }
  230.  
  231. /*************************/
  232. /*  Bit reverse counter  */
  233. /*************************/
  234.   j = 1; 
  235.   n1 = n - 1;
  236.   
  237.   for (i = 1; i <= n1; i++) 
  238.   {
  239.     if (i < j) 
  240.     {
  241.     xt = *(px+j);
  242.     *(px+j) = *(px+i); 
  243.     *(px+i) = xt;
  244.     xt = *(py+j); 
  245.     *(py+j) = *(py+i);
  246.     *(py+i) = xt;
  247.     }
  248.     
  249.     k = n / 2; 
  250.     while (k < j) 
  251.     { 
  252.     j = j - k; 
  253.     k = k / 2; 
  254.     }
  255.     j = j + k;
  256.   }
  257.  
  258. /*
  259.   for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
  260. */
  261.   
  262.   return(n);
  263. }
  264.  
  265.  
  266. /*****************************************************/
  267. /* Various timer routines.                           */
  268. /* Al Aburto, aburto@marlin.nosc.mil, 20 Dec 1992    */
  269. /*                                                   */
  270. /* t = dtime() outputs the current time in seconds.  */
  271. /* Use CAUTION as some of these routines will mess   */
  272. /* up when timing across the hour mark!!!            */
  273. /*                                                   */
  274. /* For timing I use the 'user' time whenever         */
  275. /* possible. Using 'user+sys' time is a separate     */
  276. /* issue.                                            */
  277. /*                                                   */
  278. /* Example Usage:                                    */
  279. /* [timer options added here]                        */
  280. /* main()                                            */
  281. /* {                                                 */
  282. /* double starttime,benchtime,dtime();               */
  283. /*                                                   */
  284. /* starttime = dtime();                              */ 
  285. /* [routine to time]                                 */
  286. /* benchtime = dtime() - starttime;                  */
  287. /* }                                                 */
  288. /*                                                   */
  289. /* [timer code below added here]                     */
  290. /*****************************************************/
  291.  
  292. /*********************************/
  293. /* Timer code.                   */
  294. /*********************************/
  295. /*******************/
  296. /*  Amiga dtime()  */
  297. /*******************/
  298. #ifdef Amiga
  299. #include <ctype.h>
  300. #define HZ 50
  301.  
  302. double dtime()
  303. {
  304.    double q;
  305.  
  306.    struct   tt {
  307.       long  days;
  308.       long  minutes;
  309.       long  ticks;
  310.    } tt;
  311.  
  312.    DateStamp(&tt);
  313.  
  314.    q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;
  315.  
  316.    return q;
  317. }
  318. #endif
  319.  
  320. /*****************************************************/
  321. /*  UNIX dtime(). This is the preferred UNIX timer.  */
  322. /*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
  323. /*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
  324. /*****************************************************/
  325. #ifdef UNIX
  326. #include <sys/time.h>
  327. #include <sys/resource.h>
  328.  
  329. #ifdef __hpux
  330. #include <sys/syscall.h>
  331. #define getrusage(a,b) syscall(SYS_getrusage,a,b)
  332. #endif
  333.  
  334. struct rusage rusage;
  335.  
  336. double dtime()
  337. {
  338.    double q;
  339.  
  340.    getrusage(RUSAGE_SELF,&rusage);
  341.  
  342.    q = (double)(rusage.ru_utime.tv_sec);
  343.    q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
  344.    
  345.    return q;
  346. }
  347. #endif
  348.  
  349. /***************************************************/
  350. /*  UNIX_Old dtime(). This is the old UNIX timer.  */
  351. /*  Use only if absol